Diffusion-induced instability and chaos in random oscillator networks 
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We demonstrate that diffusively coupled limit-cycle oscillators on random networks can exhibit 
various complex dynamical patterns. Reducing the system to a network analog of the complex 
Ginzburg-Landau equation, we argue that uniform oscillations can be linearly unstable with re- 
spect to spontaneous phase modulations due to diffusional coupling - the effect corresponding to 
the Benjamin-Feir instability in continuous media. Numerical investigations under this instabil- 
ity in random scale-free networks reveal a wealth of complex dynamical regimes, including partial 
- - - amplitude death, clustering, and chaos. A dynamic mean-field theory explaining different kinds of 

nonlinear dynamics is constructed. 
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_ Phase oscillators coupled through various network structures have been extensively analyzed as a prototype model 
["t I 1 of network dynamics J^, 2, 3, 4]. In most studies, complete synchronization of all oscillators has been the main focus, 
^ ■ though possibilities of more complex dynamics have also been reported 0, [1] . Coupled phase oscillators are obtained 
from general coupled limit-cycle oscillators by eliminating amplitude degrees of freedom in the weak coupling limit 0] ■ 
When the coupling is not weak, such phase reduction breaks down and much richer dynamics can be expected. 
I""!. In this paper, we analyze complex dynamics exhibited by diffusively coupled limit-cycle oscillators on random 
networks. In continuous media, sufhciently large difference in diffusion constants of oscillating components (e^. 
chemical species) can destabilize uniform oscillations and lead to diffusion- induced spatiotemporal chaotic regimes |7| , 
■ such as those experimentally observed in surface chemical reactions {S\. We argue that diffusional mobility of the 
' components can also lead to the instability and complex dynamics on networks. 
Ch Rather than treating a specific model of limit-cycle oscillators, we focus on a network version of the complex 

Ginzburg-Landau (CGL) equation derived from a general model of diffusively coupled limit-cycle oscillators near 
the supercritical Hopf bifurcation. Our linear stability analysis based on Laplacian eigenvectors of the network 
^ ' generally shows that the uniformly oscillating solution can become unstable when the analog of the Benjamin-Feir 
CS| ' (BF) condition is satisfied. Numerical simulations on random scale-free networks under this condition reveal different 
\l kinds of complex dynamical regime. To explain them, an approximate mean-field theory is constructed. 
^ [ We consider a system of diffusively coupled identical limit-cycle oscillators on random networks consisting of N 
nodes described by 
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ON : X,(t) = F(X,) + D ^ L,,Xfc. (1) 
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Here, Xj(t) represents the state of the oscillator on node j (j = f , • • • , N), F(X) specifies the intrinsic dynamics of an 
oscillator, and the last term takes into account diffusive coupling on the network, where D is a diffusion matrix and Ljk 
J-j ^ is a Laplacian matrix of the network. The network is defined by a symmetric adjacency matrix Ajk, whose components 
are I, if the nodes j and k are connected, and otherwise. The Laplacian matrix is given by Ljk = Ajk — kjSjk with 
kj = X^fcLi ^jfc representing the degree (number of connections) of node j. We assume that each oscillator has a 
stable limit-cycle solution X°(t) in absence of diffusion. A uniformly oscillating solution of the system, Xj(i) = X''(t) 

for Vj, always satisfies Eq. ([T]) because X^fcLi ^jk — holds, but diffusion may destabilize this solution. 

We assume that each oscillator is slightly above the supercritical Hopf bifurcation point and consider a situation 
where the effect of diffusion is also comparably small. Then, using the standard weakly nonlinear analysis [3], we can 
reduce Eq. (jT]) to a network version of the CGL (or Kuramoto-Tsuzuki) equation, 

N 

Wj{t) = (1 + ica)Wj - (1 + 1C2) \Wj\^ Wj + K{1 + ici) J2 LjkWk. (2) 

k=l 

Here, Wj {t) represents the complex oscillation amplitude of j-th oscillator such that Xj {t) — X^'^) oc Wj (t) exp(zti;oi)U-|- 
c.c. where X^'^^ is the unstable fixed point, ujq is the Hopf frequency, and U is the complex critical eigenvector 
of the Jacobian matrix of F(X) at X'-'^-'. Real parameters Cq, Ci, C2, and positive coupling strength K can be 
determined when F(X) and D are explicitly given. Note that if the diffusion constants of all components are equal, 
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i.e. D = DI where I is the identity matrix, we have ci = 0. Equation ([2|) has a uniformly oscillating solution, 
Wj{t) = W°{t) exp[i(co - C2)t] for Vj. 

When K is small enough, each oscillator state is always near the unperturbed limit cycle, so that Eq. ([2]) can further 
be reduced to coupled phase oscillators of the form (j)j{t) = uj — CJ2k=i ^jk sin(0j — 0fc + 7), where 4>j is the phase 
of the oscillator j, iv — cq — a is the frequen cy, C is the rescaled coupling strength, and the coupling phase shift 7 
satisfies cos 7 = (1 + C1C2)/ \/ (1 + c\){\ + c^). Recently, it has been shown that this network phase model exhibits 
coexistence of drifting and phase-locked oscillators with stationary phase gradients In the following, we focus on 
the case with stronger coupling. 

Let us analyze linear stability of the uniform solution. Plugging weakly perturbed solution Wj{t) = T/F°(i){l + 
exp[i0j(t)] into Eq. Q with Pj{t) and 9j{t) being amplitude and phase perturbations, respectively, we obtain 
the following linearized equations: 

N 
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Oj {t) = -2c2Pj +kY^ Ljk {ciPk + Ok). (3) 

To proceed, we introduce Laplacian eigenvalues A^") and eigenvectors 0'-"'' — {4''f'\ • • ■ , ^jv^) '^^ Laplacian matrix 
Ljk satisfying X^fcLi ^jk4'k'^ = A^")^^"'' for a = 1, • • • , iV. All eigenvalues are real and non-positive, and the eigen- 
vectors are mutually orthogonal. We expand the perturbations as {pj,Oj) = X]^=l(P^"^ ^'^P ('^^"^O where 
/O*-"-* and 6^"^ are expansion coefficients and A^"^ is the complex growth rate of a-th eigenmode. Then, a characteristic 

equation {A^")}^ + 2{1 - i^A^^'jA*") - 2(1 + ciC2)KA^°'^ + (1 + 4) {iCA^")}^ = is obtained for each eigenmode, 
which yields 

A^"^ = -1 + i^A'") ± y^l + 2ciC2{i^A(")} -c2{ifA(")}'. (4) 

When Re A^"^ > for some a, the a-th eigenmode is unstable. 

By expanding the upper branch A^"^ of Eq. g]) for small KA''°'\ we obtain A^"' = (1 + ciC2)if A(°) + 0({ii:A(")}^). 

Therefore, Re A^'' can be positive when the condition 1 -t- C1C2 < is satisfied (note that A^") < 0). This is the same 
as the BF condition for instability of the uniform solution of the CGL equation in continuous media Q, which also 
applies to globally coupled and non-locally coupled CGL oscillators fl^, . Note that the BF condition cannot be 
satisfied for ci = and therefore a sufficiently large difference in diffusion constants of the components is necessary. 
For the instability to actually occur, the discrete Laplacian eigenvalues should exist near the peak of the upper curve 
given by Eq. (g]) . As we already know for other coupling schemes 0, [H, , Eq- © is expected to exhibit strongly 
nonlinear behavior once the uniform solution becomes unstable. 

As an example of random networks, we use random scale-free networks of size N = 1000 and mean degrees (k) — 20 
generated by the Barabasi- Albert preferential attachment rule 0. We fix the parameters ci = — 2 and C2 = 2 (cq 
can be set to without loss of generality), and vary the coupling strength K. Numerical results shown below are for 
one particular realization of the random network, but similar behavior was observed for other network realizations as 
well. 

Figure [Ha) displays the degree kj of each node vs. the node index j, where the node indices {j} are sorted in 
decreasing order of their degrees {kj} so that inequalities fci > k2 ■ ■ ■ > k^ hold. We use this ordering as a useful 
way to visualize the complex dynamics on the network throughout our analysis. Figure [ijb) shows the Laplacian 
eigenvalues A^"^ of the same network. The eigenvalue indices {a} are also sorted in decreasing order of the eigenvalues 
such that = At^' > A^^) > • ■ • > A(^) hold. 

Figure [Ijc) plots growth rates of the perturbations Re Af^'' obtained by the linear stability analysis as functions 
of — i^TA*^"^ ai K = 0.04. The actual growth rates are distributed discretely on the curves given by Eq. ([4]). We can 
see that the growth rates on the upper branch Re A^|"'* can become positive when the coupling strength K is in an 
appropriate range, indicating that the uniform solution can undergo a diffusion-induced instability. 

To investigate nonlinear dynamics after the instability, we have performed numerical simulations of Eq. Q with 
slightly perturbed uniform solutions as initial conditions. When K is very small {K < 0.001), no oscillator deviates 
largely from the unperturbed limit-cycle orbit W'^^'>{t), so that the reduced phase model is valid. The coupling phase 
shift is given by 7 = arccos(— 3/5) ~ —2.21, which is repulsive because I7I > 7r/2 Therefore, the oscillators do not 
synchronize but rotate incoherently. When K is very large {K > 0.165), there exist no discrete growth rates on the 
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FIG. 1: (Color online) (a) Degree kj and (b) Laplacian eigenvalue A^'*^ of the scale-free network used in numerical simulations, 
(c) Linear growth rates of perturbations Re X^_^'^ plotted as functions of —KA'^"^ . Dashed straight line shows the slope —{I+C1C2) 
of the upper branch at the origin. 



positive part of the upper curve of Fig. [T]Jc), so that the uniform solution remains stable even if the BF condition is 
satisfied. 

Between these limits, we have found three characteristic steady dynamical regimes as shown in Fig. [21 where 
snapshots of the amplitude profile | Wj \ and the distribution of Wj on the complex plane are displayed for three values 
of the coupHng strength, K = 0.02, K = 0.04, and K = 0.08. 

(i) Partial amplitude death [Figs. [2l[a),(d)]. When 0.006 < K < 0.028, a group of oscillators with small node indices 
(i.e. with large degrees) stops rotation and stays near the origin of the complex plane while other oscillators are 
rotating around circular orbits incoherently, with a rather sharp but smooth transition between the two groups. 

(ii) Chaos [Figs. [2l^b),(e)]. When 0.028 < K < 0.078, the oscillators are roughly separated into three groups. In the 
first group, oscillators take approximately constant amplitudes near 0.5, which corresponds to the central cluster on 
the complex plane. Amplitudes of oscillators in the second group are strongly scattered and evolve chaotically, but 
their envelope still forms smooth curves. This group corresponds to the intermediate scattered oscillator states on 
the complex plane. The oscillators in the last group again take constant amplitudes near 1, which correspond to the 
oscillator states elongated along the unit circle on the complex plane. The largest Lyapunov exponent of the system 
is positive in this regime. 

(iii) Clustering [Figs. [2I^c),(f)]. When 0.078 < K < 0.164, phase relations among the oscillators are frozen and the 
whole system exhibits a rigid constant rotation. For relatively small values of K [K < 0.12), the oscillators with 
small degrees split into two groups with two distinct amplitudes, i.e. they exhibit a 2-cluster state. As K becomes 
larger, the two clusters gradually approach each other and, at relatively large K [K > 0.14), the two clusters merge 
to a single cluster but still with phase scattering. 

Transitions between the above dynamical regimes occur abruptly and are clearly detectable, whereas the change in 
the dynamics within each regime, e.g. transformation from 2-cluster to 1-cluster states, occurs gradually with K. 

To explain the observed dynamical patterns, we employ the mean-field approximation, valid for large random 
networks with strong diffusive mixing. It has been used in analyzing network-based epidemics spreading models , 
coupled phase oscillators [3, 5], and also network Turing patterns [ll|. A crucial point here is that we consider not 
only static but also dynamic mean fields that oscillate periodically with time. 

Introducing a complex local field hj{t) — X^feLi the diffusion term in Eq. ^ can be written as 
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FIG. 2: (Color online) (a-c): Snapshots of the complex amplitude Wj on the complex plane, (d-f): Snapshots of the amplitude 
\Wj\ vs. the node index j. Coupling strength is iC" = 0.02 for (a,d), K — 0.04 for (b,e), and K = 0.08 for (c,f). Node indices are 
sorted in decreasing order of their degrees {kj} so that inequalities fci > fc2 ■ • • > fcjv hold. Solid curves in (d-f) are predictions 
of the mean-field theory. 



LjkWk = hj{t) — kjWj. We approximate this local field as 

^ k 

h,{t)^k,H{t), H{t)^Y.T^^o(t)^ (5) 

, l^total 

where ktotai = X^^Li H{t) is a degree- weighted global mean field over the network [1, [l^ilili- Thus, we ignore 

detailed connections of the network and retain only the degrees. Equation ^ is then approximated as 

Wj(t) = (1 + ico)Wj - (1 + ic2) \Wj\^ Wj + kjK{l + ici) {H{t) - Wj} , (6) 

which describes independent CGL oscillators coupled to a global mean field H{t). The effective coupling strength of 
each oscillator to H{t) is given by kjK , and thus depends on the node degree kj. 

In Figs.lHI^a), (d), and (g), time sequences of the global mean field H{t) obtained numerically for the three cases in 
Fig. [2]are shown. H{t) almost vanishes aX K — 0.02, whereas it oscillates sinusoidally at if = 0.04 and K = 0.08. We 
can thus approximate H{t) in these regimes as H{t) = Bexp (ii7t), where B and Q, denote amplitude and frequency of 
the periodic sinusoidal oscillation, which reasonably fit the numerical data as shown in the figures. Similar sinusoidal- 
field approximation has been used in the analysis of collective dynamics of globally coupled CGL oscillators but 
degree inhomogeneity in networks essentially changes the results. Precisely speaking, in the chaotic regime. Hit) is 
only approximately sinusoidal and can be more complex, e.g. quasiperiodic for some other values of K (as also known 
in the case of global coupling [12)) but we focus on the simplest sinusoidal case here. In the clustering regime, H{t) 
is always strictly sinusoidal. 

Moving to a rotating frame by introducing W{t) — V{t) exp {iVLt), we obtain an autonomous equation for V{t) as 
V{t) — [l + i{c{) — ^l)\V ~ {1 + 102) \V\ V + (3{l + ici) [B — V). Here we dropped the index j, because all oscillators obey 
the same dynamics, and defined (3 — f3{j) = kjK, which plays the role of a bifurcation parameter. The dependence 
of the oscillator dynamics on the node index j enters only through /3. 

Figures[3]Jb) , (e), and (h) display the bifurcation diagrams of the above equation as functions of the control parameter 
/3, where the maximal and the minimal values of \W\ = \V\ are plotted using B and Q estimated numerically in 
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FIG. 3: (Color online) (a), (d), (g); Evolution of real and imaginary parts of the global mean field and fitting by B exp(jr2t). 
(b), (e), (h): Bifurcation diagrams of the sinusoidally-driven oscillator. Insets show limit cycle orbits or fixed points of V{t) at 
the parameter values indicated by broken vertical lines, (c), (f), (i): Probability density functions of \Wj\ compared with the 
mean-field approximation. Solid curves represent the maximal and minimal values of the complex amplitude \ Wj\. Parameters 
are K = 0.02, B = (a,b,c), K = 0.04, B = 0.442, Q = -1.25 (d,e,f), and K = 0.08, B = 0.532, Q, = -0.796 (g,h,i). 
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Figs. [Slja), (d), and (g). Depending on the values of B^ft, and /?, the equation exhibits a symmetric limit cycle, an 
asymmetric limit cycle, and one or two fixed points [l5|. 

Now, using the relation between the bifurcation parameter and the node degree, (3 = (3{j) = Kkj, we can map 
the bifurcation diagrams onto actual amplitude patterns in the network. The solid curves in Figs. [H^d), (e), and (f) 
are the maximal and minimal values of \Wj\ = |Vj|, which fit the envelopes of the oscillator dynamics reasonably 
well. Figure[3l^c), (f), and (i) compare the numerical probability density functions of the amplitude with these curves, 
showing good agreement. In particular, the condition for an oscillator to fall in the amplitude death state in regime (i) 
can be obtained analytically by linear stability analysis of the fixed point Vj = Q with i? = 0. This yields kj > \/K, 
which also agrees well with the numerical data. Thus, the complex network dynamics in our model can be well 
understood through the mean- field approximation p^ . 

Summarizing, we have investigated diffusion-induced instability and resulting complex dynamics exhibited by limit- 
cycle oscillators on random networks. Under the mean-field approximation, the observed inhomogeneous dynamical 
patterns can be interpreted as a mixture of various limit cycles and fixed points, which is reminiscent of the "chimera" 
states found in nonlocally coupled oscillators [1, In the present case, however, the degree inhomogeneity of the 
network essentially determines the dynamics of each oscillator. 

Dynamical systems coupled through various networks are ubiquitous structures in the real world, ranging from 
neuronal circuits in the brain to various engineering problems, such as sensor networks and power grids (see The 
fact that complex dynamical patterns can spontaneously emerge in random oscillator networks may be of fundamental 
importance in understanding the behavior and functions of such systems. 
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